Correlation Analysis of Synchronization Type and Degree in Respiratory Neural Network

Pre-Bötzinger complex (PBC) is a necessary condition for the generation of respiratory rhythm. Due to the existence of synaptic gaps, delay plays a key role in the synchronous operation of coupled neurons. In this study, the relationship between synchronization and correlation degree is established for the first time by using ISI bifurcation and correlation coefficient, and the relationship between synchronization and correlation degree is discussed under the conditions of no delay, symmetric delay, and asymmetric delay. The results show that the phase synchronization of two coupling PBCs is closely related to the weak correlation, that is, the weak phase synchronization may occur under the condition of incomplete synchronization. Moreover, the time delay and coupling strength are controlled in the modified PBC network model, which not only reveals the law of PBC firing transition but also reveals the complex synchronization behavior in the coupled chaotic neurons. Especially, when the two coupled neurons are nonidentical, the complete synchronization will disappear. These results fully reveal the dynamic behavior of the PBC neural system, which is helpful to explore the signal transmission and coding of PBC neurons and provide theoretical value for further understanding respiratory rhythm.


Introduction
Synchronous in the neuronal network, a complex population-firing pattern, is believed to play a critical role in many brain functions and many fundamental biological functions. For example, it has been implicated in sensory processing [1], the generation of sleep rhythms [2], the generation of respiratory rhythms [3,4], Parkinsonian tremor [5], epileptic seizures [6][7][8], and motor activity [9]. Synchronous arise through interactions between three network components. ese are the intrinsic properties of neurons within the network, the synaptic properties of connections between neurons, and the topology of network connectivity. Neurons with intrinsic properties can be described by a large number of mathematical models. One kind of the well-known and biologically plausible models are the Hodgkin-Huxley equations and their simplified versions, which are widely used to analyze the synchronous in the neuronal network [10][11][12]. When a large number of conducted neuron models are considered, the number of coupled differential equations can be often a problem for computer simulations. erefore, some models that are simpler but keep some of the dynamical features are considered, such as a discrete-time two-dimensional map proposed by Rulkov, and the integrate-and-fire neural models or hybrid neuron models [13,14]. Protachevicz et al. investigated how the excitatory and inhibitory connectivities from one brain area to another influence the phase angle and neuronal synchronization, in which the neuron dynamics is given by the adaptive exponential integrate-and-fire model [15]. Borges et al. also used this model to study how spiking or bursting synchronous behavior appears as a function of the coupling strength and the probability of connections [16]. Reis et al. investigated the synchronization properties of a neuronal synchronization model inspired by the connection architecture of the human cerebral cortex using the Rulkov two-dimensional discrete-time map [17]. Recently, functional neuron models were used to connect and couple neural circuits to form a reliable neural network, which could be applied to predict and regulate the collective behaviors of neural networks and neural circuits [18][19][20]. For example, the synchronization approach between thermosensitive neuron and light-dependent neuron are helpful to detect the changes in environment and heat source exactly [21]. And two piezoelectric sensing neurons (PSNs) were driven by the same external voice for detecting possible synchronization approach without any synapse coupling, which found that two identical PSNs driven by the same periodical stimuli (external forces) could reach synchronous bursting, spiking, and periodical firings, and the synchronization stability was dependent on the external forcing applied on the two PSNs in case of chaotic firing [22]. Except for functional neuron models, field coupling was also considered in the neural circuits recently. Yin et al. used Josephson junction to build a coupling channel for connecting two FitzHugh-Nagumo neural circuits; the hybrid synapse could estimate the effect of the external magnetic field by generating additive phase error between the junction [23]. Yao et al. coupled two Chua circuits via an induction coil coupling, which could benefit the realization of synchronization between two chaotic Chua systems [24].
In the early nineteenth century, a focus of what we now call neuroscience was to discover the noeud vital, that is, the site where the rhythm of breathing originates; 200 years later, the pre-Bötzinger complex (PBC) in the medulla was believed to be the kernel for breathing [25]. Network synchronized activity in the PBC in the mammalian respiratory brainstem controls the inspiratory phase of the respiratory rhythm [26]. To explore the primary kernel for respiratory rhythm in the PBC, Butera et al. presented two minimal models of the PBC and investigated the control and synchronization of coupled neurons [27,28]. Toporikova and Butera presented a two-compartment mathematical model (TB model) of isolated neurons with two independent bursting mechanisms [29]. Park and Rubin modified the TB model and proposed a single-compartment model (still known as the TB model) [30]. In order to analyze the synchronized activity in a coupled network, we simplify the TB model by considering a parameterized path in a plane, which is of the form of an ellipse with the principal axis along the intracellular calcium concentration axis and the fraction of IP3 (inositol triphosphate) channels axis, to obtain a modified TB model.
A subset of cells within the PBC are able to engage in rhythmic bursting when parameters are adjusted appropriately [31]. Using the methods of dynamical theory, the complex bursting dynamics and their transition mechanisms in the PBC neuron model have been studied widely [32][33][34][35][36]. Recently, Wang et al. considered a model of an isolated embryonic PBC neuron featuring two distinct bursting mechanisms and uncovered mechanisms underlying several different types of intrinsic bursting dynamics observed experimentally including several forms of plateau bursts, bursts involving depolarization block, and various combinations of these patterns [37]. Lü et al. found mixed bursting (MB) could also be driven by the sole action of intracellular calcium oscillations originating from the dendrite based on Park and Rubin's model for a single PBC neuron, which MB was called the dendritic mixed bursting (DMB) [38].
A single neuron exhibits complex nonlinear dynamics behavior, while a group of neurons can show more complicated patterns. In order to understand the dynamics of the neural network, the effects of various forms of coupling connection [39][40][41] and the time delay [42,43] are also researched extensively. Synchronization of inspiratory neurons in PBC has an important effect on respiratory rhythm [44]. Gaiteri and Rubin investigated how the dynamics of individual neurons (quiescent/bursting/tonic) and the betweenness centrality of neurons' positions within the network connectivity graph interact to govern network burst synchrony, by simulating heterogeneous networks of PBC neurons [45]. Ashhad and Feldman showed that the inspiratory rhythm emerges as the network reorganizes from random tonic activity toward periodic short-term synchronization [46]. Duan et al. studied the bursting dynamics of the two-coupled PBC complex neurons and explored the possible forms of dynamics that the model network could produce as well the transitions of in-and anti-phase bursting, respectively [47]. ey also observed a new type of mixed burst similar to a depolarization block bursting (DBbursting) in the model of PBC neurons and studied the types of mixed burst and their transition mechanisms by using the multi-time-scale dynamics and one-and two-parameter bifurcation analysis, as well as investigated the effects of persistent sodium conductance on the anti-phase synchronization pattern in coupled PBC neurons [48]. In this paper, using a modified TB model, the relationship between synchronization and correlation degree is established by using ISI bifurcation and correlation coefficient based on the analysis of synchronization and synchronous transition. e structure of this paper is as follows. Section 2 first introduces the modified single atrioventricular neuron model under the action of stimulation current, and then the firing interval of this neuron model is determined by using the interspike interval (ISI) bifurcation and phase plane analysis, and the chaos phenomenon in this interval is judged by combining the maximum Lyapunov exponent. At the same time, it is also suspected that the addition period of ISI sequence is related to special bifurcation, which proves that the transition from period 2 to period 3 corresponds to Hopf bifurcation. In Section 3, coupling PBC neural network is considered, and synchronization transition caused by stimulus current, coupling strength, and delay parameters is discussed in three cases: identical and no delay coupling, symmetric time-delay coupling, and asymmetric time-delay coupling. A large number of two-parameter coloring graphs are calculated, which more vividly summarizes the change rule of the synchronization process and proves the continuous dependence of delay equation on coupling strength under weak time delay. In Section 4, the most complex asymmetric nonidentical coupling is considered. In this case, it is difficult for two coupled neurons to achieve complete synchronization.
is section mainly analyzes the weak phase synchronization and calculates the similarity function 2 Computational Intelligence and Neuroscience of the membrane potential of two coupled neurons by software simulation (see Appendix for simulation environment). Finally, in Section 5, the relevant conclusions of this paper are summarized.

Single Compartment PBC Neuron Model
In 2011, Toporikova and Butera [29] established a twocompartment PBC neuron model based on experiments, and then Park and Rubin [30] modified the model to obtain a single-compartment PBC model (still known as the TB model). In this paper, an improved TB model is adopted, which is composed of two parts: the somatic subsystem and the calcium subsystem. e somatic subsystem is affected by calcium-activated nonspecific cationic (CAN) current, while the calcium subsystem is independent of the former. In this paper, we consider a parameterized path in the plane ([Ca], l), which is of the form of an ellipse with the principal axis along the [Ca] axis and l axis (see also [49] for similar ideas). e variation along the ellipse can also be regarded as the solution of the equation, as follows: [Ca] .
where [Ca] is the intracellular calcium concentration, l is the fraction of IP3(inositol triphosphate) channels in the endoplasmic reticulum membrane that have not been inactivated, [Ca] c and l c define the center of the ellipse, d is its aspect ratio, and ε is the speed with which the ellipse is traced. See Appendix for the other expression and parameter description of the function in the formula. e stimulation current I exc is an important parameter that causes PBC neurons to produce complex firing behavior. Figure 1 is the ISI bifurcation diagram of I exc . With the increase of parameters, the number of peaks in the cluster gradually increases, and the ISI sequence of neurons changes from period 1 to period 2 and then increases one by one. However, when the stimulation current increases to 11.5, the ISI sequence disappears, which means that neurons enter a resting state. As shown in Figure 2(a), I exc � 8.5, the time history diagram shows the period 3 bursting, and the phase diagram corresponding to the plane (h, V) is shown in Figure 2 It is worth noting that the threshold value of periodadding bifurcation may correspond to a special bifurcation type, and the relationship between them is discussed from the perspective of bifurcation. e bifurcation diagram of the system with respect to parameter I exc is shown in Figure 3. e solid line is the bifurcation curve of the equilibrium point, and its trend is classical cubic, in which red and black represent stable and unstable equilibrium points, respectively. e lower branch of the curve consists of a stable node and an unstable saddle point, which changes through Hopf bifurcation point (HB 2 ); then turns to the middle branch via first fold bifurcation (LF), which is an unstable equilibrium point; and finally turns to the upper branch of the equilibrium point curve through the second fold bifurcation (RF), which is basically the same as the lower branch. With the decrease of parameters, it changes from unstable to stable, and the limit cycle bifurcation occurs at HB 1 , in which green and blue represent stable and unstable periodic orbits, respectively, and the limit cycle changes from unstable to stable through saddle node on the invariant cycle (SNIC) and finally disappears into the equilibrium bifurcation curve, and the end point is homoclinic bifurcation (HC). According to this idea, the limit cycle bifurcation can also be extended in HB 2 (I exc � 8.073), and the enlarged figure is shown in the upper right corner. Although the structure type is basically the same as that at HB 1 , it is a small change, and there are pitchfork bifurcation (BP, I exc � 8.0805) and period doubling bifurcation (PD, I exc � 8.08193), which is similar to the period in Figure 1. It can be judged that the change of neuron peak has a great relationship with the type of bifurcation point, especially HB 2 . Now calculate the first Lyapunov coefficient at this point to judge the direction of Hopf bifurcation as follows: Computational Intelligence and Neuroscience C(x, y) � e equilibrium point of the original system at HB 2 is x HB � (− 51.8239, 0.003315, 0.6824, 0.1, 0.9), but because the calcium subsystem is independent, we only need to fix [Ca] at the equilibrium ([Ca] � 0.1) and discuss the somatic subsystem. e Jacobian matrix at the point is (6), and the expression of function F � (f 1 , f 2 , f 3 ) T is shown in Appendix. e matrix has a pair of conjugate complex roots λ 1,2 � ± ωi, ω � 0.0023475; the other eigenvalue is λ � − 0.869856; and the eigenvector corresponding to the positive imaginary root is Take another vector p satisfying A T p � − ωip and 〈p, q〉 � 1 and get A is the Jacobian matrix at the equilibrium point, which can be written as follows: where B(x, y) and C(x, y, z) are multiple-linear functions with component form, and the expression is where ξ � (ξ 1 , ξ 2 , ξ 3 ) is the equilibrium of the system, which can be obtained from equations (7) and (8).
It can be calculated accordingly: e first Lyapunov coefficient is erefore, the direction of Hopf bifurcation is supercritical. e detailed derivation of the first Lyapunov coefficient can be referred to reference [50].

Chaotic Coupled PBC Model
e dynamic behavior of two electrically coupled (linearly coupled) PBC neurons consists of eight differential equations including Eqs. (4) and (5) and the following equations: where i � 1, 2 stands for neurons 1 and 2. Because the two subsystems have unidirectional influence, they act on the same calcium subsystem (4, 5) when coupling two neurons. is section mainly studies the synchronization behavior of two identically coupled chaotic PBC neurons. By calculating the ISI bifurcation, correlation coefficient, and maximum synchronization difference of the two coupled neurons, it shows the rich synchronization transition modes of PBC neurons. Fixed I exc1 � I exc2 � 8.5, the correlation coefficient R and the maximum synchronization difference max (e) are defined as follows: where N is the total number of samples, V i is the membrane potential of the neuron i, V i is the corresponding average value, and ‖ · ‖ ∞ is the infinite norm.

Complete and Incomplete Synchronization.
In this section, the synchronous transition caused by coupling strength is mainly analyzed when there is no delay coupling, as shown in equation (15).
Here, g c is the coupling strength. When two neurons start to couple, a single neuron will also show different firing Computational Intelligence and Neuroscience patterns, as shown in Figure 4, which shows the relationship between the ISI of the first neuron and the coupling strength g c . When g c ∈ (− 0.5, − 0.3), neuron 1 will always period-3 bursting, which is the same as the action potential of neurons in Figure 2. With the increase of coupling strength, the size of the attractor of coupled neurons tends to increase. In this case, the coupling strength of the chaotic system destroys the attractor with three cycles in a single neuron, resulting in more complex attractors. When the two neurons are positively coupled, this situation becomes more obvious, and the peak value in the corresponding bursting becomes more irregular.
In addition, the ISI bifurcation diagram has obvious high and low parts. e higher ISI sequence corresponds to the resting state, and it will gradually decrease with the increase of coupling strength.
e coupling strength is not only an important factor affecting the firing patterns of neurons but also a key parameter affecting the synchronous transition of neural networks. It can be observed from Figures 5(a) and 5(b) that when the coupling strength is between (− 0.5, − 0.3), the correlation coefficient of membrane potential between the two neurons is always 1, and the maximum synchronization difference is always 0, which means that the two coupled neurons are completely synchronized, but with the increase of the coupling strength, they start to become incompletely synchronized. Complete synchronization and incomplete synchronization are two extreme cases, which can be directly judged by a correlation coefficient, while there are other types of synchronization under incomplete synchronization. As shown in Figure 6(a), when g c � − 0.5, the phase diagram of the coupled system on the (V 1 , V 2 ) plane coincides with the 45-degree bisector, which again indicates that the coupled neurons are in-phase synchronous state, that is, complete synchronization. With the increase of coupling strength, the degree of synchronization gradually weakens (the correlation coefficient drops sharply in Figure 5(a)). In Figure 6(b), when g c � − 0.24, the action potentials of the two coupled neurons are almost identical, but in the phase diagram on the (V 1 , V 2 ) plane, a small disturbance around the 45-degree bisector can be clearly observed, which indicates that the firing sequences of the two neurons are almost coincident but also completely overlapped, so it is called almost in-phase synchronous state.
In addition, in Figure 5(b), the maximum value of synchronization difference jumps rapidly with the increase of coupling strength to g c � − 0.3, which also indicates that complete synchronization is lost and remains basically unchanged in a certain range. When the positive coupling is achieved, the maximum synchronization difference tends to increase monotonously. If the synchronization difference is nonzero, it only means that the two neurons are not fully synchronized, and it is not certain that the coupled neurons are out of synchronization, such as phase synchronization, lag synchronization, and chaotic synchronization, which may occur. As shown in Figure 7(a), the periodic bursting of two coupled neurons occur at the same time, and the neurons generate irregular electricity, and the corresponding phase diagram is in an irregular state, and there is no fixed phase between the two firing sequences, indicating that the coupled neurons are in an asynchronous state (Figure 7(c)).
When g c � 0.4, the phase diagram is symmetrical about the 45-degree bisector, and the coupled neurons are out-ofphase synchronization. Both neurons produce almost the same bursting, but the oscillation time of action potential is inconsistent, which is one of the neuron action potentials, another in the resting state, and the time history diagram of poor synchronization also shows periodic phenomenon (Figure 7(f )). In the respiratory system, when the expiratory neurons produce exhalation, the inspiratory neurons are usually resting, and when the inspiratory neurons operate, the expiratory neurons are also in the resting state. In the respiratory system, when the expiratory neurons produce exhalation, the inspiratory neurons are usually resting, and when the inspiratory neurons operate, the expiratory neurons are also in the resting state. is phenomenon has also been found in other biological behaviors, such as the walking of limb animals, the formation of a central pattern generator (CPG), and the alternating arm swinging behavior of runners.

e Continuous Dependence on the Coupling Strength under Weak Delay.
In the real biological system, because the information transmission is carried out at a limited speed, the time delay is inevitable in the process of information transmission, which makes the dynamic system exhibit rich dynamic characteristics, such as phase locking, synchronization, and multi-stability, which shows that the time delay is of great significance in information processing based on neural movement.
In equation (15), where τ is the time delay, and i, j ∈ 1, 2 { }, and i ≠ j. is section mainly discusses the synchronization behavior caused by time delay. According to the relationship between two coupled neurons without time delay in the case of identical coupling, the synchronization effect of time delay on the coupled system is considered when the coupling strength g c � − 0.5 and g c � 0.4, respectively. As before, the infinite norm of correlation coefficient and synchronization difference is still used to judge the synchronization between coupled neurons. When g c � − 0.5, according to definitions (18) and (19), the relationship between the correlation coefficient and the maximum synchronization difference with time delay is shown in Figure 8. When the time delay is 0, it is the situation discussed above. e two coupled neurons are completely synchronized and remain in a complete synchronization state with the increase of the small range of time delay. When the delay increases to 10, the maximum synchronization difference is not 0 and basically remains unchanged, and the corresponding correlation coefficient is also very small. When the correlation coefficient is less than 0.3, it can basically be considered that there is no correlation between the two statistics, which can also be used to explain that the two coupled neurons are not completely synchronized at τ < 10. Combined with Figure 8(c), the amplitudes of the two neurons in this range are also inconsistent, thus judging that the two neurons are asynchronous. However, when the time delay is large, the correlation coefficient is close to 1, and the corresponding maximum synchronization difference fluctuates in a low amplitude state, which means that the two neurons transform between complete synchronization and approximate synchronization and are infinitely close to complete synchronization. When the τ ∈ (10, 14), the fluctuation of correlation coefficient increases and the corresponding maximum synchronization difference gradually decreases, indicating that the asynchronous state disappears and the coupled neurons reach the transitional stage of synchronization.   Computational Intelligence and Neuroscience In order to better observe the effect of the time delay on the synchronization, the release sequences of two neurons are obtained when τ 1 � 5, as shown in Figure 9(a). e release sequences of two neurons is consistent basically, and resting and active stage of bursting synchronization is achieved. However, through the difference of synchronization in Figure 9(c), it is not difficult to find that in the firing sequence of membrane potential, the firing time of two coupling neurons is different, which indicates that the two neurons are bursting synchronization with unrelated peaks, and the situations in Figures 9(b) and 9(c) are completely consistent, indicating that two coupling neurons remain chaotic when the delay is increased. Figure 9(d) is the phase diagram of the system on the (V 1 , V 2 ) plane when τ � 17, which coincides with the 45 degree, which means that the two neurons are completely synchronized.
When g c � 0.4, the ISI sequences of the two neurons are as shown in Figure 10, which consists of two parts. e higher sequence is the resting period, and the resting state increases slightly with the change of time delay. e lower sequence corresponds to the firing state of the action potential, and the peak number in the bursting changes obviously in the middle lag. On the whole, the ISI sequences of the two neurons are deviated, indicating that the membrane potential is not completely synchronized. Figure 10 shows that the system will show chaos under time delay coupling, which is also illustrated by the maximum Lyapunov exponent in Figure 11, but the chaos at this time is completely caused by time delay, instead of the transition from period to chaos in Figure 4. In this complicated situation, the correlation coefficient and the maximum potential difference are still used to judge. As shown in Figure 12, the maximum synchronization difference is kept around 65 mV, and the corresponding correlation coefficient is also small, and the fluctuation decreases with the increase of τ. It will be shown that although the two coupling neurons are not completely synchronized in this case, they can achieve weak synchronization.
When discussing identical coupling no delay, the synchronization phenomenon caused by coupling strength g c has been analyzed. When g c � − 0.5, the coupled neurons are completely synchronized, which changes from no delay to symmetric time delay. Whether there is a mutation or continuous dependence on parameters needs further analysis. In Figure 8, it is found that in the case of weak delay (τ ∈ (0, 0.02)), the coupling neurons have a high degree of synchronization, which can be understood as the continuous dependence of the equation on the parameters. e continuous dependence is mainly for ordinary differential equations, but the time delay system is not fully sure. Now, we analyze the dependence of the system on the coupling strength under weak delay.
First, the coupled system is considered as follows: where ), and f is a nonlinear differentiable function. e common method to study the synchronization problem of the coupled system is to deal with the synchronization error, which can be obtained from (20) by e system is linearized at e � 0; we have where M �  For the system used in this paper, τ 1 � τ 2 � τ ≠ 0, under weak time delay, x(τ − t) Taylor expansion is obtained as follows: Since τ is in a small range, the higher-order term in the formula can be ignored and can be obtained equation (24). Hence, Linearize at e � 0 and get where _ e � (1 + τC)(Df(u) − M)e. From the point of view of dynamics, the stability of the system (20) depends only on M, that is, C. On the other hand, the synchronization dynamics of a system with weak time delay can be transformed into the stability of the system (20), so the synchronization in this small range can be determined by the coupling strength.
e correlation coefficient and the maximum synchronization difference are observed in the two-parameter space (τ, g c ) ( Figure 13). When g c > 0, the correlation coefficient is low, and the maximum synchronization difference is too large. At this time, for any time delay, the coupled system will not achieve complete synchronization, which means that the time when the two neurons' membrane potential firing sequence is always inconsistent in the case of positive coupling, and it will not achieve bursting synchronization with the change of τ, which is obvious in negative coupling. It is mainly divided into two parts, τ < 10 and τ > 10. When τ > 10, the correlation coefficient is on the high side, and the transition from approximate synchronization to complete synchronization is gradually realized with the influence of coupling strength. When τ < 10, the correlation coefficient is low, and the maximum synchronization difference is always nonzero, which indicates that the two coupling neurons are asynchronous. In the local neighborhood of time delay and coupling strength, their similarity functions by defined (30) are mostly 0 (Figure 13(c)), indicating that neurons are completely synchronized.
is phenomenon can also be explained as the continuous dependence of the delay equation on parameters, which is completely consistent with the previous derivation.

e Relationship between Phase Synchronization and Weak Correlation under Asymmetric Time-Delay Coupling.
For different neurons, the speed of information transmission is also different, so it is necessary to discuss the asymmetric time-delay coupling system. e emergence of asymmetric time delays makes the finite-dimensional system become an infinite-dimensional system, thus inducing more complex dynamic characteristics.
where τ 1 ≠ τ 2 , i, j ∈ 1, 2 { }, and i ≠ j. e correlation coefficient of two neurons coupled with asymmetric time delays on the two-parameter plane (τ 1 , τ 2 ) is shown in Figure 14, and the correlation degree of the diagonal part is obviously higher than that of other places, which is caused by the complete synchronization caused by the strong symmetric time delay discussed in the previous section, which again conforms to the continuous dependence of the time-delay equation on parameters mentioned above, and this is applicable not only to weak time delay but also to strong time delay. However, in the nondiagonal region, the correlation coefficient is obviously low, which means that the complete synchronization of the two coupled neurons is lost from symmetric time delay to asymmetric time-delay coupling. When the correlation coefficient R � 1, it means that the coupling neurons are completely synchronized. When the correlation coefficient is high but not equal to 1, it indicates that the neurons may be from a certain degree of approximate synchronization to complete synchronization. When the correlation coefficient is high but not equal to 1, it indicates that the neurons may be from a certain degree of approximate synchronization to complete synchronization. e former case of incomplete synchronization is mainly judged by phase plane analysis, which has obvious limitations and cannot explain the characteristics in the general state.
A complete synchronization is a kind of synchronization behavior achieved by two identical chaotic systems, but the interacting chaotic systems usually are not identical in real life, which makes it difficult for coupled chaotic systems to achieve complete synchronization, but other synchronizations cannot be excluded. erefore, it is necessary to study weak synchronization behavior, and phase synchronization is one of them. In order to study phase synchronization, it is necessary to define the phase of chaotic systems. At present, there are many known methods, and the commonly used method is an analytical signal approximation and Poincaré mapping, there we use the Poincaré mapping definition phase.
where t n is the time for the trajectory to cross the Poincaré graph for the n time, if the phase difference Δϕ � |ϕ 1 (t) − ϕ 2 (t)| < const, the two coupled neurons are said to be in phase synchronization. is paper limits const ≤ 2π.
First, fixed g c � − 0.5, as shown in Figure 14, the correlation coefficient in the nondiagonal region shows a symmetrical decreasing trend, and the vertical line τ 1 � 5 is taken in the figure, and the phase difference about τ 2 is as shown in Figure 15, most of which lies above the baseline (2π), mainly corresponding to the part where R < 0.1 in Figure 14, which indicates that the coupled neurons are asynchronous, and the diagram corresponding to the phase plane is similar to that in Figures 7 and 9. Jumping near the  baseline means that neurons transform between asynchronous and phase synchronization. When the correlation coefficient is greater than 0.2 (i.e., τ 2 > 14), the two coupled neurons are always in phase synchronization. When τ 1 � 10, the coupled neurons were still within a certain range of asynchronous, but the scope was obviously smaller than τ 1 � 5, and then they are always in phase synchronization. When τ 1 increases to 17, this situation is more prominent, and it is asynchronous in a small range of time-delay changes, while the remaining ranges are all phase synchronization. e Poincaré mapping is used to analyze the period and frequency. Figure 16 shows the period and frequency of different τ 1 , and red and blue correspond to cell 1 and cell 2, respectively. Obviously, the period and frequency are inversely proportional. In Figure 16(a), when τ 2 > 14, the period and frequency of the two neurons are almost coincident, which is consistent with the phase synchronization obtained when τ 1 � 5 in Figure 15, while when τ 2 < 14, this is the same. Although the change trend is the same, the deviation is large, which means that the two neurons are asynchronous.
When τ 1 � 10, the period and frequency (Figure 16(b)) of the two coupled neurons are almost consistent with the conclusion in Figure 14, and the coupled neurons produce certain errors in some areas. e figure shows that the period and frequency do not coincide and then start to coincide after a sharp decline, and the variation amplitude is obviously reduced, which corresponds to the phase synchronization when τ 1 � 10 in Figure 14

Lag Synchronization under Nonidentical Coupling
In the previous section, the influence of coupling strength and time delay on the model under identical coupling was discussed. In this section, the synchronization phenomenon of two coupled neurons under nonidentical coupling was mainly discussed. For the coupled nonidentical chaotic system, if the coupling strength is weak, the amplitudes of the coupled neurons are irrelevant, and the chaotic system will achieve phase-locking. With the increase of the coupling strength, the amplitudes will start to build a relationship, and there will be a lag synchronization phenomenon, that is, the two states are almost the same, only the time difference exists, and one system lags behind the other. At this time, the lag synchronization can be observed through the dynamic characteristics of the amplitudes. When the coupling strength exceeds a certain critical value, the system states will remain almost the same. When the coupling strength exceeds a certain critical value, the state of the system will keep almost the same, but there will be a time delay τ, which makes x(t + τ) ≈ y(t). In order to describe the lag synchronization, a similar function is introduced as follows [51]: where x, y represents the average value, and if S(τ) is sufficiently small, it means that the two chaotic systems have achieved lag synchronization with a lag time of τ. In which a represents the average value, and if S(τ) is sufficiently small, it means that the two chaotic systems have achieved lag synchronization with a time delay τ. is section mainly studies the synchronous transition of neurons under nonidentical and asymmetric time-delay coupling. First, we consider the relationship between the firing state of coupled neurons. At that time, Figure 17 correspond to neurons 1 and 2, respectively. It can be observed from the figure that under negative coupling strength, the ISIs of two coupled neurons are almost the same, but with the increase of coupling strength, they show great differences, especially. In fact, under the condition of nonidentical coupling, the amplitude of two neurons is different, and it is also obvious in frequency. Figure 17(b) shows the relationship between the similar function S(0) and the coupling strength. It can be clearly seen that in the positive coupling, S(0) monotonically increases, while in the negative coupling, the value of S(0) is relatively small, that is, the correlation degree of the membrane potential of two neurons is relatively low, which indicates that the two coupled neurons have achieved lag synchronization. Second, the firing sequence of neurons with nonidentical coupling is significantly different, which is much more complex than that of a single neuron, which means that the stimulation current is an important factor affecting the firing pattern of neurons and also affects the synchronous behavior. As shown in Figures 17(a) and 17(b), the coloring diagram of S(0) in the parameter space (I exc1 , I exc2 ) under different coupling strengths is irregular in both cases; the correlation degree of coupling strength g c � − 0.5 is higher than that of g c � 0.4. e results show that different stimulus currents are easier to achieve lag synchronization under negative coupling. Finally, the synchronization effect of time delays on coupled neurons is discussed, and the asymmetric time delays are discussed from two aspects. On the one hand, the similar function between time delays and coupling strength g c is shown in Figure 18(c). e effect of time delays is obviously greater than that of coupling strength, and S(0) gradually increases with the increase of time delays, indicating that a single time delay cannot effectively increase the degree of synchronization, and the effect of negative coupling is more obvious than that under positive coupling. On the other hand, given the coupling strength g c � − 0.5, the relationship between similarity function and asymmetric time delays is shown in Figure 18(d), which means that the lag increases at any time and the degree of synchronization is gradually enhanced, which indicates that the time delays effectively increases the lay synchronization of coupled neurons.

Conclusion
In this paper, the synchronization of two linearly coupled PBC neurons is analyzed. By calculating the ISI, correlation coefficient, maximum synchronization difference, and similarity function of the coupled chaotic system, it is found that there is a complex synchronous transition behavior in the coupled PBC network, including asynchronization, weak synchronization, and complete synchronization. e degree of synchronization between the two coupling neurons is judged by the correlation coefficient. When the correlation coefficient is 1, it can be seen that the coupling neurons are completely synchronized. When the correlation coefficient is not 1, it is incomplete synchronous. is situation contains a variety of possible weak synchronization and requires further analysis.
Several common types of synchronization are discussed based on the effect of time delay on coupled systems in this study. First, the synchronization and transition caused by coupling strength are discussed under the condition of no delay coupling. It is found that when the coupling is negative (the coupling symbol is determined by the direction of the current), the two coupling neurons undergo fully synchronous period-3 bursting, and then they get transit to approximate synchronization. With the increase of coupling strength, the coupled neurons realize the transition from asynchronous to out-of-phase synchronization. en, symmetric time-delay coupling is discussed. Numerical simulation shows that two coupled neurons are still in complete synchronization with weak time delay. It is proved that the continuous dependence of time-delay equation on coupling strength with weak time delay combined with theory. With the increase of time delay, coupled neurons turn to an asynchronous state, and when time delay increases to a certain extent, coupled neurons begin to move repeatedly between complete synchronization and approximate synchronization (Figure 8). is synchronous transition type is more complex than the weak to the strong synchronous transition of chaotic coupled Morris-Lecar neurons [11]. In addition, on the basis of symmetric time delay, the effects of different coupling strengths are discussed, and it is found that the synchronous transition phenomenon caused by negative coupling is obviously richer than that caused by positive coupling. en, the asymmetric time-delay coupling is discussed, and the phase difference is defined by Poincaré mapping [13]. e results show that the firing pattern and synchronous transition of chaotic coupled PBC neurons are robust to large time delays, and the robustness is more obvious with the increase of time delays, that is, the asymmetric time delays always show the phase synchronization of chaotic firing in a large range. Robustness means that some performance indexes of the system remain unchanged under disturbance. At the beginning of the article, the influence of external stimulus current on a single neuron has been discussed and combined with bifurcation theory (Figure 3); it is proved that the bifurcation curve of periodic orbit generated from HB bifurcation is an important factor causing the increase of the peaks number. Finally, based on the previous related results, the nonidentical asymmetric time-delays coupling is discussed, and the delayed synchronization caused by stimulus current, time delays, and coupling strength is analyzed with the help of similarity function. e results show that the neurons are more likely to achieve lay synchronization under negative coupling and further show that the increase of time delays effectively improves the degree of synchronization of coupled neurons. In contrast, the role of stimulation current in synchronization is very complex, which needs further study. Generally speaking, with the increase of coupling strength, the synchronization transition process of two nonidentical coupled chaotic systems is from phase synchronization to lag synchronization and then to almost complete synchronization [51]. Although there are differences in the transfer process in this PBC coupling system, it also shows that the coupling strength can cause a complex synchronization process.
Physiological experiments have shown that neuronal synchronization is closely related to many clinical diseases [5,6]. erefore, it is necessary to discuss the synchronous transition of neurons. is study shows that even if it is a simple linear coupling, by controlling the coupling strength and time delay, all kinds of synchronization can be achieved effectively.
For convenience, make f 1 (V, n, h) � − I L + I K + I Na + I NaP + I CAN + I exc C m , ere are many parameters in this model including Eq. (4), Eq. (5), and Eqs. (15)(16)(17), and the values used in this paper are listed in Table 1.
f is the hill function related to calcium concentration: e description of the somatic subsystem can be referred to references [29,30], and the derivation and explanation of the imposed path about the calcium subsystem are shown in reference [49]. e numerical software in this paper is mainly XPPAUT and MATLAB, and the fourth-order Runge-Kutta algorithm is used with a step size of 0.1.

Data Availability
e simulation data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest. Table 1: Parameter values used for the functions in the eightdimensional model equations (4)(5) and (15)(16)(17).

Parameter
Value C m 21 μF g Na 28 nS g K 11.2 nS g L 2.3 nS g NaP 2 n g CAN 0.7 nS V Na 50 mV